www.gusucode.com > Non-photorealistic Camera工具箱源码matlab程序 > Non-photorealistic Camera/PIE.m

    function im_out = PIE( im_target,lambda,c)
%PIE function: blends the source image with the target one based on the
%boundary given as a BW mask using Poisson Image Editing (PIE)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Initialization

if size(im_target,1)~= size(lambda,1) || ...
        size(im_target,2)~= size(lambda,2)
    error('Please, use images with the same size');
end

%if the mask is not grayscale, convert it
im_mask=ones(size(lambda));
im_mask(1:3,:)=0; im_mask(end-3:end,:)=0;
im_mask(:,1:3)=0; im_mask(:,end-3:end)=0;



%convert source and target images to double for more precise computations
im_target=double(im_target);

lambda=double(lambda);

%if we are working with true color, let m=3 otherwise, m=1
if c==0 %true color images
    c=3; %for the next for loop
    
else
    c=1;
    if size(im_target,3)>1
        im_target=rgb2gray(im_target);
    end
    
end

%initially, output image = target image
im_out=im_target;

%normalize the mask image to assure that unknown pixels = 1
im_mask=mat2gray(im_mask);

%convert it to logical to ignore any fractions (soft masks)
im_mask=im_mask==1;

%find the number of unknown pixels based on the mask
n=size(find(im_mask==1),1);

%create look up table
map=zeros(size(im_mask));

%loop through the mask image to initialize the look up table for mapping
counter=0;
for x=1:size(map,1)
    for y=1:size(map,2)
        if im_mask(x,y)==1 %is it unknow pixel?
            counter=counter+1;
            map(x,y)=counter;  %map from (x,y) to the corresponding pixel
            %in the 1D vector
        end
    end
end


for i=1:c %for each color channel
    
    % loop through the mask image again to:
    %1- initialize the coefficient matrix
    %2- initialize the B vector
    
    %if the method is seamless cloning; so, intially put B= (-) laplacian of
    %im_source,
    %otherwise (mixing gradients), B= (-) max(laplacian of im_source, laplacian
    %of im_target)
    
    %create the coefficient matrix A
    
    %At most, there are 5 coefficients per row according to eq (3)
    %in the report
    coeff_num=5;
    
    %create the sparse matrix to save memory
    A=spalloc(n,n,n*coeff_num);
    
    %create the right hand side of the linear system of equations (AX=B)
    B=zeros(n,1);
    
    %create the gradient mask for the first derivative
    grad_mask_x=[-1 1];
    grad_mask_y=[-1;1];
    
    %get the first derivative of the target image
    g_x_target=conv2(im_target(:,:,i),grad_mask_x, 'same');
    g_y_target=conv2(im_target(:,:,i),grad_mask_y, 'same');
    

    g_x_final=g_x_target.*(lambda);
    g_y_final=g_y_target.*(lambda);
    
    %get the final laplacian of the combination between the source and
    %target images lap=second deriv of x + second deriv of y
    lap=conv2(g_x_final,grad_mask_x, 'same');
    lap=lap+conv2(g_y_final,grad_mask_y, 'same');
    

counter=0;
for x=1:size(map,1)
    for y=1:size(map,2)
        if im_mask(x,y)==1
            counter=counter+1;
            A(counter,counter)=4; %the diagonal represent the current pixel
            
            %check the boundary
            if im_mask(x-1,y)==0 %known left pixel
                B(counter)=im_target(x-1,y,i); %add it to B
            else %unknown boundary
                A(counter,map(x-1,y))=-1; %set its coefficient to -1
            end
            if im_mask(x+1,y)==0 %known right pixel
                B(counter)=B(counter)+im_target(x+1,y,i); %add it to B
            else %unknown boundary
                A(counter,map(x+1,y))=-1; %set its coefficient to -1
            end
            if im_mask(x,y-1)==0 %known bottom pixel
                B(counter)=B(counter)+im_target(x,y-1,i); %add it to B
            else %unknown boundary
                A(counter,map(x,y-1))=-1; %set its coefficient to -1
            end
            if im_mask(x,y+1)==0 %known top pixel
                B(counter)=B(counter)+im_target(x,y+1,i); %add it to B
            else %unknown boundary
                A(counter,map(x,y+1))=-1; %set its coefficient to -1
            end
            
            %update the B vector with the laplacian value
            
            B(counter)=B(counter)-lap(x,y);
            
        end
    end
end

%solve the linear system of equation
X=A\B;


%reshape X to restore the output image

%     counter=0;
%     for x=1:size(map,1)
%         for y=1:size(map,2)
%             if im_mask(x,y)==1
%                 counter=counter+1;
%                 im_out(x,y,i)=X(counter);
%             end
%         end
%     end


for counter=1:length(X)
    [index_x,index_y]=find(map==counter);
    im_out(index_x,index_y,i)=X(counter);
    
end


%release all
clear A B X lap_source lap_target g_mag_source g_mag_target
end
im_out=real(im_out)/255;


end